      subroutine area_bc(i1, i2, j1, j2, ks, iwid, jwid, kwid, updown, x, y, z&
         , ax, ay, az) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: i1 
      integer , intent(in) :: i2 
      integer , intent(in) :: j1 
      integer , intent(in) :: j2 
      integer , intent(in) :: ks 
      integer , intent(in) :: iwid 
      integer , intent(in) :: jwid 
      integer , intent(in) :: kwid 
      real(double) , intent(in) :: updown 
      real(double) , intent(in) :: x(*) 
      real(double) , intent(in) :: y(*) 
      real(double) , intent(in) :: z(*) 
      real(double) , intent(out) :: ax(*) 
      real(double) , intent(out) :: ay(*) 
      real(double) , intent(out) :: az(*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: j, i, ijk 
      real(double) :: x1, x2, y1, y2, z1, z2 
!-----------------------------------------------
!
      do j = j1, j2 
         do i = i1, i2 
!
            ijk = (ks - 1)*kwid + (j - 1)*jwid + (i - 1)*iwid + 1 
!
            x1 = 0.5*(x(ijk+iwid)+x(ijk+iwid+jwid)-x(ijk)-x(ijk+jwid)) 
!
            x2 = 0.5*(x(ijk+jwid)+x(ijk+iwid+jwid)-x(ijk)-x(ijk+iwid)) 
!
            y1 = 0.5*(y(ijk+iwid)+y(ijk+iwid+jwid)-y(ijk)-y(ijk+jwid)) 
!
            y2 = 0.5*(y(ijk+jwid)+y(ijk+iwid+jwid)-y(ijk)-y(ijk+iwid)) 
!
            z1 = 0.5*(z(ijk+iwid)+z(ijk+iwid+jwid)-z(ijk)-z(ijk+jwid)) 
!
            z2 = 0.5*(z(ijk+jwid)+z(ijk+iwid+jwid)-z(ijk)-z(ijk+iwid)) 
!
            ax(ijk) = (y1*z2 - y2*z1)*updown 
            ay(ijk) = (z1*x2 - z2*x1)*updown 
            az(ijk) = (x1*y2 - x2*y1)*updown 
!
         end do 
      end do 
!
      return  
      end subroutine area_bc 


      subroutine chek_surf(i1, i2, j1, j2, ksurf, iwid, jwid, kwid, x, y, z, px&
         , py, pz, up, vp, wp, ax, ay, az, dt, nlist, ijktmp3) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: i1 
      integer , intent(in) :: i2 
      integer , intent(in) :: j1 
      integer , intent(in) :: j2 
      integer , intent(in) :: ksurf 
      integer , intent(in) :: iwid 
      integer , intent(in) :: jwid 
      integer , intent(in) :: kwid 
      integer , intent(out) :: nlist 
      real(double) , intent(in) :: px 
      real(double) , intent(in) :: py 
      real(double) , intent(in) :: pz 
      real(double) , intent(in) :: up 
      real(double) , intent(in) :: vp 
      real(double) , intent(in) :: wp 
      real(double) , intent(in) :: dt 
      integer , intent(out) :: ijktmp3(*) 
      real(double) , intent(in) :: x(*) 
      real(double) , intent(in) :: y(*) 
      real(double) , intent(in) :: z(*) 
      real(double) , intent(in) :: ax(*) 
      real(double) , intent(in) :: ay(*) 
      real(double) , intent(in) :: az(*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: i, j, ijk 
      real(double) :: testout, testin 
!-----------------------------------------------
!
!     a routine to make a list of those polygons the particle orbit
!     may have intersected
!
!
      nlist = 0 
!
      do i = i1, i2 
         do j = j1, j2 
!
            ijk = (ksurf - 1)*kwid + (j - 1)*jwid + (i - 1)*iwid + 1 
!
            testout = (px - x(ijk))*ax(ijk) + (py - y(ijk))*ay(ijk) + (pz - z(&
               ijk))*az(ijk) 
!
            testin = (up*ax(ijk)+vp*ay(ijk)+wp*az(ijk))*dt 
!
!test    ********************************************
!
            if (testout<0. .or. testout>=testin) cycle  
!
            nlist = nlist + 1 
            ijktmp3(nlist) = ijk - kwid 
!
         end do 
      end do 
!
      return  
      end subroutine chek_surf 


      subroutine closest_vertex(i1, i2, j1, j2, k1, k2, iwid, jwid, kwid, x, y&
         , z, xp, yp, zp, imin, jmin, kmin, ijkmin) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: i1 
      integer , intent(in) :: i2 
      integer , intent(in) :: j1 
      integer , intent(in) :: j2 
      integer , intent(in) :: k1 
      integer , intent(in) :: k2 
      integer , intent(in) :: iwid 
      integer , intent(in) :: jwid 
      integer , intent(in) :: kwid 
      integer , intent(out) :: imin 
      integer , intent(out) :: jmin 
      integer , intent(out) :: kmin 
      integer , intent(out) :: ijkmin 
      real(double) , intent(in) :: xp 
      real(double) , intent(in) :: yp 
      real(double) , intent(in) :: zp 
      real(double) , intent(in) :: x(*) 
      real(double) , intent(in) :: y(*) 
      real(double) , intent(in) :: z(*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: k, j, i, ijk 
      real(double) :: dmin, d 
!-----------------------------------------------
!
!
      dmin = 1.E10 
!
      do k = k1, k2 
         do j = j1, j2 
            do i = i1, i2 
!
               ijk = (k - 1)*kwid + (j - 1)*jwid + (i - 1)*iwid + 1 
!
               d = sqrt((xp - x(ijk))**2+(yp-y(ijk))**2+(zp-z(ijk))**2) 
!
               if (d >= dmin) cycle  
!
               dmin = d 
               imin = i 
               jmin = j 
               kmin = k 
               ijkmin = ijk 
!
            end do 
         end do 
      end do 
!
      return  
      end subroutine closest_vertex 


      subroutine look_evrywher(i1, i2, j1, j2, k1, k2, iwid, jwid, kwid, eps, &
         succes, x, y, z, xp, yp, zp, pxi, peta, pzta) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: i1 
      integer , intent(in) :: i2 
      integer , intent(in) :: j1 
      integer , intent(in) :: j2 
      integer , intent(in) :: k1 
      integer , intent(in) :: k2 
      integer  :: iwid 
      integer  :: jwid 
      integer  :: kwid 
      real(double)  :: eps 
      real(double)  :: xp 
      real(double)  :: yp 
      real(double)  :: zp 
      real(double) , intent(out) :: pxi 
      real(double) , intent(out) :: peta 
      real(double) , intent(out) :: pzta 
      logical , intent(out) :: succes 
      real(double)  :: x(*) 
      real(double)  :: y(*) 
      real(double)  :: z(*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: k, j, i, ijk, iter, itmax 
      real(double) :: nu, tsi, eta 
!-----------------------------------------------
!
!
      succes = .FALSE. 
!
      l1: do k = k1, k2 
         do j = j1, j2 
            do i = i1, i2 
!
               ijk = (k - 1)*kwid + (j - 1)*jwid + (i - 1)*iwid + 1 
!
               call map3d (ijk, iwid, jwid, kwid, eps, iter, itmax, x, y, z, xp&
                  , yp, zp, tsi, eta, nu) 
!
               if (.not.(tsi*(1. - tsi)>=0. .and. eta*(1.-eta)>=0. .and. nu*(1.&
                  -nu)>=0.)) cycle  
!
               succes = .TRUE. 
!
               pxi = i + tsi 
               peta = j + eta 
               pzta = k + nu 
!
               exit  l1 
!
            end do 
         end do 
      end do l1 
!
      return  
      end subroutine look_evrywher 


      subroutine parbctor(ncellsp, ijkctmp, iphead, npart, iwid, jwid, kwid, &
         ibar, jbar, kbar, cdlt, sdlt, x, y, z, uptilde, vptilde, wptilde, px, &
         py, pz, up, vp, wp, pxi, peta, pzta) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: ncellsp 
      integer , intent(in) :: npart 
      integer  :: iwid 
      integer  :: jwid 
      integer  :: kwid 
      integer , intent(in) :: ibar 
      integer , intent(in) :: jbar 
      integer , intent(in) :: kbar 
      real(double) , intent(in) :: cdlt 
      real(double) , intent(in) :: sdlt 
      integer , intent(in) :: ijkctmp(*) 
      integer , intent(in) :: iphead(*) 
      real(double)  :: x(*) 
      real(double)  :: y(*) 
      real(double)  :: z(*) 
      real(double) , intent(inout) :: uptilde(*) 
      real(double) , intent(inout) :: vptilde(*) 
      real(double)  :: wptilde(*) 
      real(double) , intent(inout) :: px(0:npart) 
      real(double) , intent(inout) :: py(0:npart) 
      real(double)  :: pz(0:npart) 
      real(double) , intent(inout) :: up(0:npart) 
      real(double) , intent(inout) :: vp(0:npart) 
      real(double)  :: wp(0:npart) 
      real(double) , intent(inout) :: pxi(0:npart) 
      real(double) , intent(inout) :: peta(0:npart) 
      real(double) , intent(in) :: pzta(0:npart) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: ibp1, jbp1, kbp1, n, ijk, np, inew, jnew, knew 
      real(double) :: nu, eps, wspx, wspy, wsup, wsvp 
      logical :: succes 
!-----------------------------------------------
!
!
      ibp1 = ibar + 1 
      jbp1 = jbar + 1 
      kbp1 = kbar + 1 
!
      eps = 1.E-4 
!
!
      do n = 1, ncellsp 
!
         ijk = ijkctmp(n) 
         np = iphead(ijk) 
!
         inew = int(pxi(np)) 
         jnew = int(peta(np)) 
         knew = int(pzta(np)) 
!
  380    continue 
!
         inew = int(pxi(np)) 
         if ((ibp1 - inew)*(inew - 2) < 0) then 
            if (inew < 2) pxi(np) = pxi(np) + float(ibar) 
            if (inew > ibp1) pxi(np) = pxi(np) - float(ibar) 
            go to 380 
         endif 
!
!
  390    continue 
!
         jnew = int(peta(np)) 
!
         if ((jbp1 - jnew)*(jnew - 2) >= 0) cycle  
!
         if (jnew < 2) then 
!
            peta(np) = peta(np) + float(jbar) 
!
            wspx = cdlt*px(np) - sdlt*py(np) 
            wspy = sdlt*px(np) + cdlt*py(np) 
            px(np) = wspx 
            py(np) = wspy 
!
            wsup = cdlt*up(np) - sdlt*vp(np) 
            wsvp = sdlt*up(np) + cdlt*vp(np) 
            up(np) = wsup 
            vp(np) = wsvp 
!
            wsup = cdlt*uptilde(ijk) - sdlt*vptilde(ijk) 
            wsvp = sdlt*uptilde(ijk) + cdlt*vptilde(ijk) 
            uptilde(ijk) = wsup 
            vptilde(ijk) = wsvp 
!
         endif 
!
         if (jnew > jbp1) then 
!
            peta(np) = peta(np) - float(jbar) 
!
            wspx = cdlt*px(np) + sdlt*py(np) 
            wspy = (-sdlt*px(np)) + cdlt*py(np) 
            px(np) = wspx 
            py(np) = wspy 
!
            wsup = cdlt*up(np) + sdlt*vp(np) 
            wsvp = (-sdlt*up(np)) + cdlt*vp(np) 
            up(np) = wsup 
            vp(np) = wsvp 
!
            wsup = cdlt*uptilde(ijk) + sdlt*vptilde(ijk) 
            wsvp = (-sdlt*uptilde(ijk)) + cdlt*vptilde(ijk) 
            uptilde(ijk) = wsup 
            vptilde(ijk) = wsvp 
!
         endif 
!
         go to 390 
!
      end do 
!
      return  
      end subroutine parbctor 


      subroutine parbccyl(ncellsp, ijkctmp, iphead, npart, iwid, jwid, kwid, &
         ibar, jbar, kbar, dz, x, y, z, px, py, pz, up, vp, wp, pxi, peta, pzta&
         ) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: ncellsp 
      integer , intent(in) :: npart 
      integer  :: iwid 
      integer  :: jwid 
      integer  :: kwid 
      integer , intent(in) :: ibar 
      integer , intent(in) :: jbar 
      integer , intent(in) :: kbar 
      real(double) , intent(in) :: dz 
      integer , intent(in) :: ijkctmp(*) 
      integer , intent(in) :: iphead(*) 
      real(double)  :: x(*) 
      real(double)  :: y(*) 
      real(double)  :: z(*) 
      real(double)  :: px(0:npart) 
      real(double) , intent(inout) :: py(0:npart) 
      real(double)  :: pz(0:npart) 
      real(double)  :: up(0:npart) 
      real(double)  :: vp(0:npart) 
      real(double)  :: wp(0:npart) 
      real(double) , intent(inout) :: pxi(0:npart) 
      real(double) , intent(inout) :: peta(0:npart) 
      real(double) , intent(in) :: pzta(0:npart) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: ibp1, jbp1, kbp1, n, ijk, np, inew, jnew, knew 
      real(double) :: nu, eps 
      logical :: succes 
!-----------------------------------------------
!
!
      ibp1 = ibar + 1 
      jbp1 = jbar + 1 
      kbp1 = kbar + 1 
!
      eps = 1.E-4 
!
      do n = 1, ncellsp 
!
         ijk = ijkctmp(n) 
         np = iphead(ijk) 
!
         inew = int(pxi(np)) 
         jnew = int(peta(np)) 
         knew = int(pzta(np)) 
!
  380    continue 
!
         inew = int(pxi(np)) 
         if ((ibp1 - inew)*(inew - 2) < 0) then 
            if (inew < 2) pxi(np) = pxi(np) + float(ibar) 
            if (inew > ibp1) pxi(np) = pxi(np) - float(ibar) 
            go to 380 
         endif 
!
!
  390    continue 
!
         jnew = int(peta(np)) 
!
         if ((jbp1 - jnew)*(jnew - 2) >= 0) cycle  
!
         if (jnew < 2) then 
!
            peta(np) = peta(np) + float(jbar) 
!
            py(np) = py(np) + dz 
!
         endif 
!
         if (jnew > jbp1) then 
!
            peta(np) = peta(np) - float(jbar) 
!
            py(np) = py(np) - dz 
!
         endif 
!
         go to 390 
!
      end do 
!
      return  
      end subroutine parbccyl 


      subroutine parlocat(ncellsp, ijkctmp, iphead, iwid, jwid, kwid, nsampl, &
         vrms, t_wall, ico, ijktmp2, ijktmp3, ijktmp4, rmaj, dz, alfa, ax, ay, &
         az, dt, itdim, npart, ibar, jbar, kbar, mgeom, cdlt, sdlt, plost_pos, &
         plost_neg, wate, x, y, z, bxv, byv, bzv, xptilde, yptilde, zptilde, &
         uptilde, vptilde, wptilde, tsix, tsiy, tsiz, etax, etay, etaz, nux, &
         nuy, nuz, link, iplost, qpar, px, py, pz, up, vp, wp, pxi, peta, pzta&
         , cartesian, killer, wl, wr, wb, wt, we, wf, dx, dy, dzdummy, xl, xr, &
         yb, yt, ze, zf, nu_len, nu_comm) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: ncellsp 
      integer  :: iwid 
      integer  :: jwid 
      integer  :: kwid 
      integer  :: nsampl 
      integer  :: itdim 
      integer  :: npart 
      integer  :: ibar 
      integer  :: jbar 
      integer  :: kbar 
      integer  :: mgeom 
      integer , intent(inout) :: iplost 
      integer  :: wl 
      integer  :: wr 
      integer  :: wb 
      integer  :: wt 
      integer  :: we 
      integer  :: wf 
      real(double) , intent(in) :: rmaj 
      real(double)  :: dz 
      real(double)  :: dt 
      real(double)  :: cdlt 
      real(double)  :: sdlt 
      real(double) , intent(in) :: plost_pos 
      real(double) , intent(in) :: plost_neg 
      real(double)  :: dx 
      real(double)  :: dy 
      real(double)  :: dzdummy 
      real(double)  :: xl 
      real(double)  :: xr 
      real(double)  :: yb 
      real(double)  :: yt 
      real(double)  :: ze 
      real(double)  :: zf 
      real(double)  :: nu_len 
      logical , intent(in) :: cartesian 
      logical  :: nu_comm 
      integer  :: ijkctmp(*) 
      integer  :: iphead(*) 
      integer  :: ico(0:npart) 
      integer , intent(inout) :: ijktmp2(*) 
      integer  :: ijktmp3(*) 
      integer  :: ijktmp4(*) 
      integer , intent(inout) :: link(0:npart) 
      integer  :: killer(*) 
      real(double)  :: vrms(*) 
      real(double)  :: t_wall(*) 
      real(double)  :: alfa(*) 
      real(double)  :: ax(*) 
      real(double)  :: ay(*) 
      real(double)  :: az(*) 
      real(double)  :: wate(itdim,*) 
      real(double)  :: x(*) 
      real(double)  :: y(*) 
      real(double)  :: z(*) 
      real(double)  :: bxv(*) 
      real(double)  :: byv(*) 
      real(double)  :: bzv(*) 
      real(double)  :: xptilde(*) 
      real(double)  :: yptilde(*) 
      real(double)  :: zptilde(*) 
      real(double)  :: uptilde(*) 
      real(double)  :: vptilde(*) 
      real(double)  :: wptilde(*) 
      real(double)  :: tsix(*) 
      real(double)  :: tsiy(*) 
      real(double)  :: tsiz(*) 
      real(double)  :: etax(*) 
      real(double)  :: etay(*) 
      real(double)  :: etaz(*) 
      real(double)  :: nux(*) 
      real(double)  :: nuy(*) 
      real(double)  :: nuz(*) 
      real(double) , intent(in) :: qpar(0:npart) 
      real(double)  :: px(0:npart) 
      real(double)  :: py(0:npart) 
      real(double)  :: pz(0:npart) 
      real(double)  :: up(0:npart) 
      real(double)  :: vp(0:npart) 
      real(double)  :: wp(0:npart) 
      real(double)  :: pxi(0:npart) 
      real(double)  :: peta(0:npart) 
      real(double)  :: pzta(0:npart) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: ibp1, jbp1, kbp1, n, nchek, itry, nchek_nxt, ijk, np, iold, &
         jold, kold, ijkold, iter, itmax, ichange, jchange, kchange, inew, jnew&
         , knew, ijknew, nout, iaxis, jaxis, nout_nxt, nlist, nl, ifail, ll, &
         jsurf, ijksurf 
      real(double) :: nu, eps, tsi, eta, epsn, updown, xi, dice, plost, cross, &
         phisurf, parphinp, parphin, zta 
      logical :: succes 
!-----------------------------------------------
!   E x t e r n a l   F u n c t i o n s
!-----------------------------------------------
      real(double) , external :: ranf 
!-----------------------------------------------
!
!
!
      if (cartesian) then 
         call fence_box (ncellsp, ijkctmp, iphead, px, py, pz, up, vp, wp, xl, &
            xr, yb, yt, ze, zf, wl, wr, wb, wt, we, wf, killer, npart) 
         call map_box (ncellsp, ijkctmp, iphead, px, py, pz, pxi, peta, pzta, &
            ibar, jbar, kbar, dx, dy, dz, npart, nu_len, nu_comm) 
         return  
      endif 
 
!
      ibp1 = ibar + 1 
      jbp1 = jbar + 1 
      kbp1 = kbar + 1 
!
!
      ijktmp2(:ncellsp) = ijkctmp(:ncellsp) 
!
!
      nchek = ncellsp 
!
!
!
!     ******************************************************************
!
!     calculate new natural coordinates
!
!     ******************************************************************
      itry = 0 
!
!
    1 continue 
      itry = itry + 1 
!
!
      nchek_nxt = 0 
!
      do n = 1, nchek 
!
         ijk = ijktmp2(n) 
         np = iphead(ijk) 
         if (np == 0) cycle  
!
         iold = int(pxi(np)) 
         iold = min(iold,ibp1) 
         iold = max(2,iold) 
!
         jold = int(peta(np)) 
         jold = min(jold,jbp1) 
         jold = max(2,jold) 
!
         kold = int(pzta(np)) 
         kold = min(kold,kbp1) 
         kold = max(2,kold) 
!
!
         ijkold = (kold - 1)*kwid + (jold - 1)*jwid + (iold - 1)*iwid + 1 
!
         eps = 1.E-4 
!
         call map3d (ijkold, iwid, jwid, kwid, eps, iter, itmax, x, y, z, px(np&
            ), py(np), pz(np), tsi, eta, nu) 
!
!      epsn=1.e-7
         epsn = 0.0 
         if (iter < itmax) then 
!
            ichange = 0.5*ifix(sign(1.,tsi - 1. - epsn) + sign(1.,tsi + epsn)) 
!
            jchange = 0.5*ifix(sign(1.,eta - 1. - epsn) + sign(1.,eta + epsn)) 
!
            kchange = 0.5*ifix(sign(1.,nu - 1. - epsn) + sign(1.,nu + epsn)) 
!
         else 
!
            kchange = sign(1.,2.*(ranf() - 0.5)) 
            ichange = sign(1.,2.*(ranf() - 0.5)) 
            jchange = sign(1.,2.*(ranf() - 0.5)) 
!
         endif 
         inew = iold + ichange 
         jnew = jold + jchange 
         knew = kold + kchange 
!
         ijknew = (knew - 1)*kwid + (jnew - 1)*jwid + (inew - 1)*iwid + 1 
!
         if (ijknew /= ijkold) then 
!
            nchek_nxt = nchek_nxt + 1 
            ijktmp2(nchek_nxt) = ijktmp2(n) 
!
!
            pxi(np) = inew + 0.5 
            peta(np) = jnew + 0.5 
            pzta(np) = knew + 0.5 
!
!
         else 
!
            pxi(np) = float(iold) + tsi 
            peta(np) = float(jold) + eta 
            pzta(np) = float(kold) + nu 
!
         endif 
!
!
      end do 
!
      nchek = nchek_nxt 
!
!
!      check for particles outside the domain
!
      nout = 0 
!
      do n = 1, nchek 
         ijk = ijktmp2(n) 
         np = iphead(ijk) 
!
         inew = int(pxi(np)) 
         jnew = int(peta(np)) 
         knew = int(pzta(np)) 
!
         if (.not.(inew<2 .or. inew>ibp1 .or. jnew<2 .or. jnew>jbp1 .or. knew<2&
             .or. knew>kbp1)) cycle  
!
         nout = nout + 1 
         ijktmp3(nout) = ijk 
!
      end do 
!
!     check for axis crossings
!
      do n = 1, nout 
!
         ijk = ijktmp3(n) 
         np = iphead(ijk) 
!
         inew = int(pxi(np)) 
         jnew = int(peta(np)) 
         knew = int(pzta(np)) 
!
         if (.not.(jnew>=2 .and. jnew<=jbp1 .and. knew<2)) cycle  
!
         iaxis = inew 
         jaxis = jnew 
!
!     the particle has crossed the axis
!
         call map3d_axis (2, ibp1, jnew, jnew, iaxis, jaxis, iwid, jwid, kwid, &
            eps, x, y, z, px(np), py(np), pz(np), tsi, eta, nu, succes) 
!
         if (.not.succes) then 
!
            call map3d_axis (2, ibp1, 2, jbp1, iaxis, jaxis, iwid, jwid, kwid, &
               eps, x, y, z, px(np), py(np), pz(np), tsi, eta, nu, succes) 
!
         endif 
         if (succes) then 
            pxi(np) = iaxis + tsi 
            peta(np) = jaxis + eta 
            pzta(np) = 2. + nu 
         else 
            pzta(np) = 2. 
         endif 
!
      end do 
!
!
      do n = 1, nout 
!
         ijk = ijktmp3(n) 
         np = iphead(ijk) 
!
 3810    continue 
!
         inew = int(pxi(np)) 
!
         if (inew>=2 .and. inew<=ibp1) cycle  
!
         if (inew < 2) pxi(np) = pxi(np) + float(ibar) 
         if (inew > ibp1) pxi(np) = pxi(np) - float(ibar) 
!
         go to 3810 
!
      end do 
!
!     nout is the number of particles for which intersections
!     have still not been found
!
!
!
      if (nout > 0) then 
!
         nout_nxt = 0 
!
!     calculate directed areas over k=kbp2 surface
!
         updown = 1. 
         call area_bc (2, ibp1, 2, jbp1, kbp1 + 1, iwid, jwid, kwid, updown, x&
            , y, z, ax, ay, az) 
!
!     one particle at a time, make list of surface elements
!     which the particle may have passed through
!
         do n = 1, nout 
            ijk = ijktmp3(n) 
            np = iphead(ijk) 
!
            succes = .FALSE. 
!
            knew = int(pzta(np)) 
            if (knew > kbp1) then 
!
!
               call chek_surf (2, ibp1, 2, jbp1, kbp1 + 1, iwid, jwid, kwid, x&
                  , y, z, px(np), py(np), pz(np), uptilde(ijk), vptilde(ijk), &
                  wptilde(ijk), ax, ay, az, dt, nlist, ijktmp4) 
!
!      if(nlist.gt.0) write(*,*) 'chek_surf, nlist=',nlist
!     calculate intersection with boundary in candidate cells
!
               knew = kbp1 
!
               do nl = 1, nlist 
                  ijknew = ijktmp4(nl) 
!
                  call map3d_surf (ijk, ijknew, iwid, jwid, kwid, eps, ifail, &
                     iphead, itdim, x, y, z, px(np), py(np), pz(np), uptilde(&
                     ijk), vptilde(ijk), wptilde(ijk), xi, eta, alfa(ijk)) 
!
!      write(*,*) 'np=',np,'ijknew=',ijknew,'ifail=',ifail
!
                  if (ifail == 1) cycle  
                  alfa(ijk) = alfa(ijk)/dt 
!
                  if (.not.(xi*(1. - xi)>=0. .and. eta*(1.-eta)>=0. .and. alfa(&
                     ijk)*(1.-alfa(ijk))>=0.)) cycle  
!
!     particle has been found
!
                  succes = .TRUE. 
!
                  knew = (ijknew - 1)/kwid + 1 
                  jnew = (ijknew - (knew - 1)*kwid - 1)/jwid + 1 
                  inew = (ijknew - (knew - 1)*kwid - (jnew - 1)*jwid - 1)/iwid&
                      + 1 
!
                  pxi(np) = xi + inew 
                  peta(np) = eta + jnew 
!
                  pzta(np) = kbp1 + 1 - 1.E-6 
!
                  ll = 1 
!
                  dice = ranf() 
!
                  if (qpar(np) > 0.0) then 
                     plost = plost_pos 
                  else 
                     plost = plost_neg 
                  endif 
!
                  if (dice >= plost) then 
!
                     call parrefl (ll, ijktmp3(n), iphead, npart, itdim, iwid, &
                        jwid, kwid, x, y, z, wate, bxv, byv, bzv, px, py, pz, &
                        up, vp, wp, pxi, peta, pzta) 
!     apply diffuse reflection conditions to the particles
!
!      call parrefl_diffuse(ll,ijktmp3(n),iphead,npart,itdim,
!     &     nsampl,vrms,T_wall,ico,
!     &     iwid,jwid,kwid,x,y,z,wate,bxv,byv,bzv,
!     &     px,py,pz,up,vp,wp,pxi,peta,pzta)
!     now we will deliberately lose the suckers after working so
!      hard to find them
!
                  else 
!c
                     iphead(ijk) = link(np) 
                     link(np) = iplost 
                     iplost = np 
!c
                  endif 
!
!      write(*,*) 'particle',np,'added to iplost'
!
!      write(*,*) 'particle found and reflected, np=', np
!
                  exit  
!
               end do 
!
            endif 
!
            if (succes) cycle  
!
            nout_nxt = nout_nxt + 1 
            ijktmp3(nout_nxt) = ijk 
!
         end do 
!
         nout = nout_nxt 
      endif 
!
!
!     check j=jbp2 surface for crossings
!
      if (nout > 0) then 
         nout_nxt = 0 
!
         updown = 1. 
!
         call area_bc (2, kbp1, 2, ibp1, jbp1 + 1, kwid, iwid, jwid, updown, x&
            , y, z, ax, ay, az) 
!
!
         do n = 1, nout 
!
            ijk = ijktmp3(n) 
            np = iphead(ijk) 
            jsurf = jbp1 + 1 
!     ******************************************************************
!     DIAGNOSTIC
            cross = -1. 
            ijksurf = 5*kwid + jbp1*jwid + 5*iwid + 1 
            phisurf = atan(y(ijksurf)/x(ijksurf)) 
            parphinp = atan(py(np)/(px(np)+1.E-6)) 
            parphin = atan((py(np)-vp(np)*dt)/(px(np)-up(np)*dt+1.E-6)) 
            if (parphinp<=phisurf .or. parphin>=phisurf) cycle  
!
!      if nlist>0, call periodic bc's
!
            call chek_surf (2, kbp1, 2, ibp1, jsurf, kwid, iwid, jwid, x, y, z&
               , px(np), py(np), pz(np), uptilde(ijk), vptilde(ijk), wptilde(&
               ijk), ax, ay, az, dt, nlist, ijktmp4) 
!
            succes = .FALSE. 
!
            if (nlist > 0) then 
!
               do nl = 1, nlist 
!
                  ijknew = ijktmp4(nl) 
!
                  call map3d_surf (ijk, ijknew, kwid, iwid, jwid, eps, ifail, &
                     iphead, itdim, x, y, z, px(np), py(np), pz(np), uptilde(&
                     ijk), vptilde(ijk), wptilde(ijk), zta, xi, alfa(ijk)) 
!
                  if (ifail /= 0) cycle  
!
                  alfa(ijk) = alfa(ijk)/dt 
!      write(*,*) 'loop 610, ifail=',ifail,'ijknew,np=',ijknew,np
!
                  if (.not.(xi*(1. - xi)>=0. .and. alfa(ijk)*(1.-alfa(ijk))>=0.&
                      .and. zta*(1.-zta)>=0.)) cycle  
!
                  succes = .TRUE. 
!
!      write(*,*) '610 loop, xi,zta,alfa=', xi,zta,alfa(ijk)
                  peta(np) = jsurf 
                  if (rmaj > 0.0) then 
!
                     call parbctor (1, ijktmp3(n), iphead, npart, iwid, jwid, &
                        kwid, ibar, jbar, kbar, cdlt, sdlt, x, y, z, uptilde, &
                        vptilde, wptilde, px, py, pz, up, vp, wp, pxi, peta, &
                        pzta) 
!
                  else 
!
                     call parbccyl (1, ijktmp3(n), iphead, npart, iwid, jwid, &
                        kwid, ibar, jbar, kbar, dz, x, y, z, px, py, pz, up, vp&
                        , wp, pxi, peta, pzta) 
!
                  endif 
                  exit  
!
               end do 
            endif 
!
            if (succes) cycle  
!
            nout_nxt = nout_nxt + 1 
            ijktmp3(nout_nxt) = ijk 
!
         end do 
!
         nout = nout_nxt 
      endif 
!
      if (nout > 0) then 
         nout_nxt = 0 
!
         updown = -1. 
         call area_bc (2, kbp1, 2, ibp1, 2, kwid, iwid, jwid, updown, x, y, z, &
            ax, ay, az) 
!
!
         do n = 1, nout 
!
            ijk = ijktmp3(n) 
            np = iphead(ijk) 
            jsurf = 2 
!
!      if nlist>0, call periodic bc's
!
            call chek_surf (2, kbp1, 2, ibp1, jsurf, kwid, iwid, jwid, x, y, z&
               , px(np), py(np), pz(np), uptilde(ijk), vptilde(ijk), wptilde(&
               ijk), ax, ay, az, dt, nlist, ijktmp4) 
!
            succes = .FALSE. 
!
            if (nlist > 0) then 
!
!     ******************************************************************
!     DIAGNOSTIC
               cross = -1. 
               ijksurf = 5*kwid + jwid + 5*iwid + 1 
               phisurf = atan(y(ijksurf)/x(ijksurf)) 
               parphinp = atan(py(np)/(px(np)+1.E-6)) 
               parphin = atan((py(np)-vp(np)*dt)/(px(np)-up(np)*dt+1.E-6)) 
!      write(*,*) 'particle np=',np,'has crossed j=2'
               if (parphinp<phisurf .and. parphin>phisurf) cross = 1. 
               do nl = 1, nlist 
!
                  ijknew = ijktmp4(nl) 
!
                  call map3d_surf (ijk, ijknew, kwid, iwid, jwid, eps, ifail, &
                     iphead, itdim, x, y, z, px(np), py(np), pz(np), uptilde(&
                     ijk), vptilde(ijk), wptilde(ijk), zta, xi, alfa(ijk)) 
!
                  if (ifail /= 0) cycle  
                  alfa(ijk) = alfa(ijk)/dt 
!      write(*,*) 'loop 710, ifail=',ifail,'ijknew,np=',ijknew,np
!
                  if (.not.(xi*(1. - xi)>=0. .and. alfa(ijk)*(1.-alfa(ijk))>=0.&
                      .and. zta*(1.-zta)>=0.)) cycle  
!
                  succes = .TRUE. 
!
!      write(*,*) '710 loop, xi,zta,alfa=', xi,zta,alfa(ijk)
                  peta(np) = jsurf - 1 
                  if (rmaj > 0.0) then 
!
                     call parbctor (1, ijktmp3(n), iphead, npart, iwid, jwid, &
                        kwid, ibar, jbar, kbar, cdlt, sdlt, x, y, z, uptilde, &
                        vptilde, wptilde, px, py, pz, up, vp, wp, pxi, peta, &
                        pzta) 
!
                  else 
!
                     call parbccyl (1, ijktmp3(n), iphead, npart, iwid, jwid, &
                        kwid, ibar, jbar, kbar, dz, x, y, z, px, py, pz, up, vp&
                        , wp, pxi, peta, pzta) 
!
                  endif 
                  exit  
!
               end do 
            endif 
!
            if (succes) cycle  
!
            nout_nxt = nout_nxt + 1 
            ijktmp3(nout_nxt) = ijk 
!
         end do 
!
         nout = nout_nxt 
!
      endif 
!
!     **********************************************
      if (nchek>0 .and. itry<5) go to 1 
!
      nchek_nxt = 0 
!
      return  
      end subroutine parlocat 


!
      subroutine trilinp(ncells, ijkcell, itdim, iwid, jwid, kwid, iphead, pxi&
         , peta, pzta, wate, bxv, byv, bzv, bxpn, bypn, bzpn) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: ncells 
      integer , intent(in) :: itdim 
      integer , intent(in) :: iwid 
      integer , intent(in) :: jwid 
      integer , intent(in) :: kwid 
      integer , intent(in) :: ijkcell(*) 
      integer , intent(in) :: iphead(*) 
      real(double) , intent(in) :: pxi(0:*) 
      real(double) , intent(in) :: peta(0:*) 
      real(double) , intent(in) :: pzta(0:*) 
      real(double) , intent(in) :: wate(itdim,*) 
      real(double) , intent(in) :: bxv(*) 
      real(double) , intent(in) :: byv(*) 
      real(double) , intent(in) :: bzv(*) 
      real(double) , intent(out) :: bxpn(*) 
      real(double) , intent(out) :: bypn(*) 
      real(double) , intent(out) :: bzpn(*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, ijk, np, inew, jnew, knew, ijknew 
!-----------------------------------------------
!
!
!dir$ ivdep
      do n = 1, ncells 
!
         ijk = ijkcell(n) 
         np = iphead(ijk) 
!
         inew = int(pxi(np)) 
         jnew = int(peta(np)) 
         knew = int(pzta(np)) 
!
         ijknew = (knew - 1)*kwid + (jnew - 1)*jwid + (inew - 1)*iwid + 1 
!
!
!
         bxpn(ijk) = wate(ijk,1)*bxv(ijknew+iwid) + (wate(ijk,2)*bxv(ijknew+&
            iwid+jwid)+(wate(ijk,3)*bxv(ijknew+jwid)+(wate(ijk,4)*bxv(ijknew)+(&
            wate(ijk,5)*bxv(ijknew+iwid+kwid)+(wate(ijk,6)*bxv(ijknew+iwid+jwid&
            +kwid)+(wate(ijk,7)*bxv(ijknew+jwid+kwid)+wate(ijk,8)*bxv(ijknew+&
            kwid))))))) 
!
         bypn(ijk) = wate(ijk,1)*byv(ijknew+iwid) + (wate(ijk,2)*byv(ijknew+&
            iwid+jwid)+(wate(ijk,3)*byv(ijknew+jwid)+(wate(ijk,4)*byv(ijknew)+(&
            wate(ijk,5)*byv(ijknew+iwid+kwid)+(wate(ijk,6)*byv(ijknew+iwid+jwid&
            +kwid)+(wate(ijk,7)*byv(ijknew+jwid+kwid)+wate(ijk,8)*byv(ijknew+&
            kwid))))))) 
!
         bzpn(ijk) = wate(ijk,1)*bzv(ijknew+iwid) + (wate(ijk,2)*bzv(ijknew+&
            iwid+jwid)+(wate(ijk,3)*bzv(ijknew+jwid)+(wate(ijk,4)*bzv(ijknew)+(&
            wate(ijk,5)*bzv(ijknew+iwid+kwid)+(wate(ijk,6)*bzv(ijknew+iwid+jwid&
            +kwid)+(wate(ijk,7)*bzv(ijknew+jwid+kwid)+wate(ijk,8)*bzv(ijknew+&
            kwid))))))) 
!
      end do 
!
      return  
      end subroutine trilinp 
